---
title: "R Notebook"
output: html_notebook
---
```{r}
library(tidyverse)
```

# Simulation under Gaussian outcomes: Used for computing limit experiment and estimating critical values
```{r}
stopping_time_Gaussian = function(eta, Delta, n) {
  flag = FALSE
  rho = 0
  i = 0
  while(flag == FALSE) {
    
    if ((i %% 2) == 0){
       mu_a = Delta 
    } else {
      mu_a = 0
    }
    rho = rho + (mu_a/n) + (1/sqrt(n))*rnorm(1, 0, 1) 
    i = i + 1
    
    if (abs(rho) >= 0.536357/eta){
      stopping_time = (i/n) 
      flag = TRUE
    } 
    
  }
  return(stopping_time)
}
```

# Finite n simulation using Uniformly distributed errors
```{r}
stopping_time = function(eta, Delta, n) {
  flag = FALSE
  rho = 0
  i = 0
  while(flag == FALSE) {
    
    if ((i %% 2) == 0){
       mu_a = Delta
       r = 1
    } else {
      mu_a = 0
      r = -1
    }
    rho = rho + (mu_a/n) + (1/sqrt(n))*sqrt(3)*r*runif(1, -1, 1)
    i = i + 1
    
    if (abs(rho) >= 0.536357/eta){
      stopping_time = (i/n) 
      flag = TRUE
    } 
    
  }
  return(stopping_time)
}

## sanity check
n = 15000
time_draws = replicate(10000, expr = stopping_time(1, 0, n))
critical_value = quantile(time_draws, 0.05)
critical_value
```

# Computing critical values hrough Gaussian approximation
```{r}
M = 15000
n_seq = seq(200, 2600, by = 200)
size = rep(0, length(n_seq))
critical_values = rep(0, length(n_seq))

for (i in 1:length(n_seq)){
  n = n_seq[i]
  time_draws = replicate(M, expr = stopping_time(1, 0, n))
  critical_values[i] = quantile(replicate(M, expr = stopping_time_Gaussian(1, 0, n)),
                                0.05)
  size[i] = mean((time_draws <= critical_values[i]))
  print(critical_values[i])
  print(size[i])
}
```

# Computation of size when using two-sample test
```{r}
M = 10000
gamma_star = 0.536357
n_seq = seq(200, 2600, by = 200)
size_two_sample = rep(0, length(n_seq))

for (i in 1:length(n_seq)){
  n = n_seq[i]
  time_draws = replicate(M, expr = stopping_time(1, 0, n))
  two_sample_test = gamma_star*sqrt(1/time_draws)
  size_two_sample[i] = mean(two_sample_test >= 1.96)
}
```

```{r}
size.df = data.frame(n_seq = n_seq,
                     size = size,
                     size_two_sample = size_two_sample) %>%
  pivot_longer(-n_seq,
               names_to = "type",
               values_to = "size") %>%
  mutate(type = if_else(type == "size", "Proposed test", "Two-sample test")) %>%
  mutate(type = factor(type, levels = c("Proposed test", "Two-sample test")))
size.df
```

```{r}
ggplot(size.df) +
  geom_point(aes(x = n_seq,
                 y = size,
                 color = type),
             size = 2.2) +
  geom_line(aes(x = n_seq,
                 y = size,
                 color = type),
            linetype = "longdash",
            size = 1) +
  labs(x = "Sample size (n)", y = "Size", color = "Type of test") +
  geom_hline(yintercept = 0.05, 
             #linetype = "dashed",
             size = 0.25,
             color = "blue") +
  theme(#aspect.ratio = 4/3,
        text = element_text(size=13),
        axis.text = element_text(size=13),
        axis.title = element_text(size=14),
        legend.position = c(0.87, 0.15),
        #legend.title = "Type of test"
        #plot.margin=grid::unit(c(0,0,0,0), "mm")
        ) +
  ylim(0, 0.1) 


#setwd("/Users/akarun/Library/CloudStorage/Dropbox/Optimal sequential tests/Figures")
ggsave("Size.png")
```


# Computation and plot of power envelopes
```{r}
#Power envelopes
M = 5000
delta_seq = seq(1, 35, by = 1)
power_seq = rep(0, length(delta_seq))

power = function(eta, n) {
critical_value = quantile(replicate(2*M, expr = stopping_time_Gaussian(1, 0, n)),
                          0.05)
  for (i in 1:length(delta_seq)) {
    delta = delta_seq[i]
    time_draws = replicate(M, expr = stopping_time(1, delta, n))
    power_seq[i] = mean((time_draws <= critical_value))
  }
  return(power_seq)
}


power_200 = power(1, 200)
power_500 = power(1, 500)
power_1000 = power(1, 1000)
power_2500 = power(1, 2500)
```

```{r}
power_limit = function(eta, n) {
critical_value = quantile(replicate(M, expr = stopping_time_Gaussian(1, 0, n)),
                          0.05)  
  for (i in 1:length(delta_seq)) {
    delta = delta_seq[i]
    time_draws = replicate(10000, expr = stopping_time_Gaussian(1, delta, n))
    power_seq[i] = mean((time_draws <= critical_value))
  }
  return(power_seq)
}

power_limit = power_limit(1, 20000)

power.df = data.frame(delta_seq, power_200, power_500, power_1000, power_2500, power_limit)
names(power.df) = c("delta_seq", "200", "500", "1000", "2500", "Power envelope")
```

```{r}
power.df1 = power.df %>% 
  pivot_longer(-delta_seq,
               names_to = "n",
               values_to = "Power") %>%
  mutate(n = factor(n, levels = c("200", "500", "1000", "2500", "Power envelope")))

ggplot(power.df1, aes(x = delta_seq,
                      y = Power,
                      color = n)) +
  geom_point() +
  geom_line() +
  xlab((expression(paste("Scaled treatment effect (", eta, ")")))) +
  ylab("Power") +
  theme(#aspect.ratio = 4/3,
        text = element_text(size=13),
        axis.text = element_text(size=12),
        axis.title = element_text(size=13),
        legend.position = c(0.85, 0.3)
        ) +
  scale_color_brewer(palette = "RdGy")

#setwd("/Users/akarun/Library/CloudStorage/Dropbox/Optimal sequential tests/Figures")
ggsave("Power envelope.png")
```